/*
 * File:   densitySP.cpp
 * Author: gian
 *
 * Created on December 12, 2012, 4:50 PM
 */

#include <math.h>
#include <iostream>
#include "densitySP.hpp"

densitySP::densitySP(int lsizex,
                     int lsizey,
                     int lsizez,
                     int nx,
                     int ny,
                     int nz) {
  mThreadNumber    = nx * ny * nz;
  mThreadArguments = new measure_thread_args[mThreadNumber];

  int deltax = lsizex / nx;
  int deltay = lsizey / ny;
  int deltaz = lsizez / nz;
  int i      = 0;

  for (int ix = 0; ix < nx; ix++) {
    for (int iy = 0; iy < ny; iy++) {
      for (int iz = 0; iz < nz; iz++) {
        mThreadArguments[i].String  = new std::vector<double>;
        mThreadArguments[i].Density = new std::vector<double>;
        mThreadArguments[i].xmin    = ix * deltax;
        mThreadArguments[i].xmax    = (ix + 1) * deltax - 1;
        mThreadArguments[i].ymin    = iy * deltay;
        mThreadArguments[i].ymax    = (iy + 1) * deltay - 1;
        mThreadArguments[i].zmin    = iz * deltaz;
        mThreadArguments[i].zmax    = (iz + 1) * deltaz - 1;

        i++;
      }
    }
  }
}

densitySP::~densitySP() {
  for (int i = 0; i < mThreadNumber; i++) {
    delete mThreadArguments[i].String;
    delete mThreadArguments[i].Density;
  }

  delete[] mThreadArguments;
}

void densitySP::placca(nodeSP *              nd1,
                       nodeSP *              nd2,
                       nodeSP *              nd3,
                       nodeSP *              nd4,
                       std::vector<double> & String,
                       std::vector<double> & Density) {
  int n = String.size();

  for (int i = 0; i < n; i++) {
    String[i] = 0.0;
  }

  latticeSP::measure_link(nd1, nd2, String);
  latticeSP::measure_link(nd2, nd3, String);
  latticeSP::measure_link(nd3, nd4, String);
  latticeSP::measure_link(nd4, nd1, String);

  for (int i = 0; i < n; i++) {
    if (String[i] > (M_PI / 2)) {
      Density[i] += 1.0;
    } else if (String[i] < (-M_PI / 2)) {
      Density[i] += 1.0;
    }
  }
}

void densitySP::operator ()(latticeSP & l) {
  mDensity.resize(l.get_N());
  mString.resize(l.get_N());

  for (int i = 0; i < l.get_N(); i++) {
    mDensity[i] = 0.0;
  }

  for (int i = 0; i < mThreadNumber; i++) {
    mThreadArguments[i].l = &l;

    pthread_create(mThreads + i, NULL, densitySP::measure_thread, mThreadArguments + i);
  }

  for (int i = 0; i < mThreadNumber; i++) {
    pthread_join(mThreads[i], NULL);

    for (int j = 0; j < l.get_N(); j++) {
      mDensity[j] += (*(mThreadArguments[i].Density))[j];
    }
  }

  for (int i = 0; i < l.get_N(); i++) {
    mDensity[i] = mDensity[i] / (3 * l.get_lattice_size());
  }
}

void * densitySP::measure_thread(void * args) {
  measure_thread_args * p = (measure_thread_args *) args;

  p -> Density -> resize(p -> l -> get_N());
  p -> String -> resize(p -> l -> get_N());

  for (int i = 0; i < p -> l -> get_N(); i++) {
    (*(p -> Density))[i] = 0.0;
  }

  for (int z = p -> zmin; z <= p -> zmax; z++) {
    for (int y = p -> ymin; y <= p -> ymax; y++) {
      for (int x = p -> xmin; x <= p -> xmax; x++) {
        nodeSP * nd   = p -> l -> operator ()(x, y, z);
        nodeSP * ndx  = p -> l -> operator ()(x + 1, y, z);
        nodeSP * ndxy = p -> l -> operator ()(x + 1, y + 1, z);
        nodeSP * ndy  = p -> l -> operator ()(x, y + 1, z);
        nodeSP * ndyz = p -> l -> operator ()(x, y + 1, z + 1);
        nodeSP * ndz  = p -> l -> operator ()(x, y, z + 1);
        nodeSP * ndxz = p -> l -> operator ()(x + 1, y, z + 1);

        placca(nd, ndx, ndxy, ndy, *(p -> String), *(p -> Density));
        placca(nd, ndy, ndyz, ndz, *(p -> String), *(p -> Density));
        placca(nd, ndz, ndxz, ndx, *(p -> String), *(p -> Density));
      }
    }
  }

  pthread_exit(NULL);
}

double densitySP::color_density(int n) {
  if ((n >= 0) && (n < mDensity.size())) {
    return mDensity[n];
  } else {
    return 0.0;
  }
}

double densitySP::total_density() {
  double total_density = 0.0;

  for (int i = 0; i < mDensity.size(); i++) {
    total_density += mDensity[i];
  }

  return total_density;
}


//~ Formatted by Jindent --- http://www.jindent.com
